Numerical study of sediment scour at meander flume outlet of boxed culvert diversion work

Background Sediment scour at downstream of hydraulic structures is one of the main reasons threatening its stability. Several soil properties and initial input data have been studied to investigate its influence on scour hole geometry by both physical and numerical models. However, parameters of resistance affecting sedimentation and erosion phenomena have not been carried out in the literature. Besides, the auxiliary-like wing walls prevalently used in many real applications have been rarely addressed for their effect on morphological change. Results In this study, a 3D Computational Fluid Dynamics model is utilized to calibrate the hydraulic characteristics of steady flow going through the culvert by comparison with experimental data, showing good agreement between water depth, velocity, and pressure profiles at the bottom of the boxed culvert. The results show that a grid cell of 0.015 m gave minimum NRMSE and MAE values in test cases. Another approach is numerical testing sediment scour at a meander flume outlet with a variety of roughness/d50 ratio (cs) and diversion wall types. The findings include the following: cs = 2.5 indicates the close agreement between the numerical and analytical results of maximum scour depth after the culvert; the influence of four types of wing wall on the geometrical deformation including erosion at the concave bank and deposition at the convex bank of the meander flume outlet; and two short headwalls represent the best solution that accounts for small morphological changes. Conclusions The influence of the roughness parameter of soil material and headwall types on sediment scour at the meander exit channel of hydraulic structure can be estimated by the numerical model.


Introduction
In dam projects, river diversion is considered the first step in the dam construction [1]. The issue is studied for three types of conveyance structures: i) tunnels and culverts, ii) open channels, and iii) overspill structures in the permanent works. Diversion works are temporary structures with short lifetimes, so their stability is often not a primary concern. However, water released into a river from tunnels or culverts may cause scouring of the riverbed, downstream of the cofferdam [1]. Therefore, the optimal design of the diversion system guarantees workshop safety and prevents possible damages resulting from conveying water during construction. The knowledge of anticipated local scours geometry has been the main concern of engineers or researchers for years because it is a significant criterion for the proper design of sluice outlet foundation [2][3][4][5]. Hence, predicting local scours after water conveyance structures such as spillways, outlet works, etc., has been widely studied to discover adequate protection solutions for the construction.
Numerous researchers have investigated the sediment scour problem. Physical models, artificial intelligence (AI) approaches, and numerical models have been efficient tools in this issue. Physical models are usually used to build empirical equations to calculate maximum scour depth and scour hole geometry [2,[6][7][8][9][10][11]. However, physical models also exposed several limitations including time-consuming and costly. Especially, it is not flexible or easy to change the dimension or to install auxiliary work as well as the initial conditions, and boundary conditions during experimenting. Besides, the narrow range of physical conditions causes limitations when applying these empirical equations in case studies. On the other hand, some Artificial Intelligence (AI) approaches have been recently applied to predict the maximum scour depth in hydraulic structures [12,13]. A study by [14] showed that AI approaches were a suitable platform to reach the scour depth prediction with a permissible level of accuracy rather than empirical equations. However, to the best of our knowledge, AI approaches have not been proposed to predict the location of the maximum scour depth and other scour hole geometries [12]. In addition to the above approaches, numerical models are known to be more flexible, thus avoiding the limitations of physical models [11,15,16]. The numerical result also gives a comprehensive representation of hydraulic features of the whole computational domain. A number of well-known Computational Fluid Dynamics (CFD) models including OpenFoam, Ansys Fluent, Flow 3D, etc., have been widely utilized in this field [8,[17][18][19]. These numerical models based on the coupling of the Volume Fluid Method and Navier Stokes equations have played important roles in simulating sediment scour issues due to the help of state-of-the-art 3D CFD models [10,20]. The robustness and effectiveness of commercial Flow 3D software are proven in practical applications [21,22]. In addition, Flow-3D was a useful software to simulate the local scour since it gives a significant association between the simulation and experiment [23,24]. However, roughness/d 50 (c s ) standing for the resistance of bed material in Flow 3D is not well studied in those studies although this parameter affects significantly the numerical solution of the local scour [8]. Additionally, headwalls at culvert outlets have been used prevalently in controlling outflow to direct in one main direction. The construction of headwalls against the flow causes a difference in velocity distribution that will generate different scouring mechanisms at downstream flumes. Previous works often study the effect of this device on flow regimes in exit channels (e.g., water depth and velocity distributions). However, there are only a few physical studies that addressed headwall design issues and local scour phenomena affected by headwalls [2,4,25]. The effect of wing walls and the effect of bed resistance parameters are underexplored in the previous numerical studies. Therefore, in this study, we intend to analyze the effect of the headwall on the geometry of scour hole after culvert as well as the erosion and deposition at the meander exit channel of diversion work using this numerical model. Besides, the hydraulic characteristics of meandering flumes are usually much more complex than those of straight flumes due to the occurrence of Coriolis force and whirlpool flow [3]. Therefore, three objectives are investigated: validate the presented model (Flow 3D) by comparing numerical results with observed data and analyze mesh sensitivity with three mesh resolutions; study the most reasonable parameter c s in simulating maximum scour depth (d s max ); discover the effects of wing wall types on the morphological pattern at the outlet of culvert and exit meander channel.
Some novel investigations can be listed: The numerical model Flow 3D is an efficient tool to simulate flow over culvert and bed deformation on exit meander channel. c s = 2.5 is a reasonable parameter in evaluating d s max . Two long headwalls cause the deepest scour hole while two short walls create the smallest sand removal volume. The meandering of the exit channel induces erosion at the concave bank and deposition at the convex bank. of the physical model of Dong Nai diversion work in the hydraulic laboratory of Thuyloi University, Vietnam with a ratio of 1/60 to examine its conveyance capacity. The diversion work including a two-box culvert, approach, and exit channels conveyed water from the upstream to the downstream of coffer dams. The intake was 117 cm in length, 11 cm in height, and 8 cm in width. Prismatic upstream and exit flumes were trapezoidal cross-sections with a 19 cm-bed width and 1.5 cm of side slope. Two wing walls at the entrance of the culvert were 68 cm in length, 13 cm in height, and 2 cm in width. Conversely, the long wall at the culvert installed on the concave bank was 91 cm in length, 13 cm in height, and 2 cm in width, while the short wall on the convex side had dimensions that were 67 cm in length, 5.0 cm in height, and 2 cm in width. The exit channel had a meander with a radius R = 125 cm [26].

Physical model
The purpose of the physical model was to investigate the conveyance capacity of the boxed culvert in both dry and rainy seasons. In the flood season, water was conducted by both temporal intake and spillways, while in the dry season, water was drained out only by the culvert. In this study, three operating conditions incorporated with the dry season were selected (Table 1).
In fact, the physical model did not carry out sediment issues on the exit channel because of the poorly measuring instruments in the laboratory. The measurement and observation were only done with pressure heads along the box tubes of the culvert, water level and velocity at upstream, and downstream sections, respectively [26]. In the second phase of this project, a mathematical model was required to assess hydraulic performance including velocity, vorticity, water depth profiles at cross sections, and sediment transport in the whole conveyance structure.

Governing equations
The commercial software package Flow 3D was built to solve Navier-Stokes equations by the Volume of Fluid (VOF) method, which is based on the conservation of two masses and momentums. Recently, this model has been increasingly used because it can be applied to simulate many sophisticated hydrodynamic problems.
Flow 3D provides two turbulent models to describe turbulent flow: Reynolds Averaged Navier Stokes (RANS) and Large Eddy Simulation (LES). The renormalization group (RNG), which is one of the RANS types, was selected due to its capacity to simulate flow over hydraulic works [27]. Moreover, the deformation of bed geometry can be demonstrated by the sediment scour module. Ten non-cohesive soil species, which are different in terms of mass density, main grain size, entrainment and transport coefficients, and critical shear stress can be added to the model as multilayers of bed information. This model can simulate the sediment transport process, which includes settling, packing, advection, bedload transport, entrainment, and depositions for each species, using the Flow 3D user manual [27].

2.2.1.Continuity equation.
The continuity equation is established by mass conservation. In general, bedload and suspended transport are used to describe the movement of sand particles in fluid flow. In the mathematical model, the bed boundary can be considered as a packed one if the local scour occurred at that place. The morphology of the packed boundary is estimated based on the conservation of mass. This process includes bedload transport, absorption, and deposition. The suspended sediment is estimated by sediment concentration and is considered a constraint at each computational cell. For each soil type, this term is estimated by the following continuity equation: where V f is volume fraction; ρ is fluid density; u, v, and w are velocity components in the x, y, and z directions, respectively; and A x , A y , and A z are the area fractions.

Momentum equations.
Three momentum equations in the x, y, and z directions are as follows: G x , G y , and G z are the body accelerations, and f x , f y , and f z are the viscous accelerations.

Sediment scour model.
The sediment transport process is often described in two phases, namely bedload transport and suspended transport. Bedload transport illustrates the motion of soil particles, such as rolling, hopping, and sliding along the packed bed surface due to the shear stress. Bedload transport means the movement of sand particles along the bed channel, regardless of whether some of them become suspended movement. The empirical formulas estimating bedload transport applied in Flow 3D were Meyer-Peter Muller, Nielsen, or Van Rijin [28,29].
The critical Shields parameter θ cr is used to define the critical bed shear stress τ cr , at which sediment movement begins for both entrainment and bedload transport, which is applied to the horizontal bed.
The Soulsby-Whitehouse equation is used to estimate the critical shear stress as follows: where ν is the kinematic viscosity of the fluid, ρ s is the soil mass density, and ρ is the fluid mass density.
The suspended sediment concentration is calculated by solving the following equation: where C s is the suspended sediment mass concentration, which is defined as the sediment mass per volume of fluid-sediment mixture; K is the diffusivity; and u s is the sediment velocity.
The computational domain in Flow 3D demonstrated the physical model of the Dong Nai diversion culvert, which was divided into two blocks (Fig 2): block 1 including the intake, upstream flume, and apron, was solid; and block 2 was a tail flume, which is defined by a packed sediment boundary. The upper boundary of block 1 and the lower boundary of block 2 had specific pressures corresponding to the water elevation of the three operating conditions, as shown in Table 1. All y-axis boundaries of blocks 1 and 2 were walls. Two objectives were set up as follows: • Simulate hydraulic features such as flow regime, water depths, and velocities at the upstream and downstream of the intake, pressure head at the gauges on the two-box tubes of the culvert, so that the selected numerical model could be validated and the mesh sensitivity could be estimated. Only 100 seconds were required for computational time to maintain a steady flow.
• Numerically estimate the impact of alternative types of wing walls on the dimensions of scour holes after culvert, four models (A, B, C, and D) of wing walls were introduced, as seen in Table 2 and Fig 3. Model A was installed in the physical model. The other cases were scenarios established using a mathematical method. Additionally, to approach the steady flow in the whole domain of the active sediment to scour module, at least 900 s-1000 s of computational time should be set up. This point will be explained in section 3.2. Based on the criterion of non-cohesive materials after hydraulic work of TCVN 9160 combined with the magnitude of velocity at the outlet of this case study, uniform soil properties were selected (i.e., entrainment coefficient of 0.005, and mass density δ s = 2650 kg/m 3 ) [30]. The mean grain size (d 50 ) chosen in this study was 3.6 mm.

Validation and mesh sensitivity analysis
To validate the proposed numerical model, flow regime, water depth, and magnitude of velocity at the inlet of intake pressure heads along each culvert's tube corresponding with the three test cases were predicted and compared with observed data. Besides, to verify mesh sensitivity, three mesh resolutions were selected to simulate hydraulic characteristics for the A model, namely 0.01 m, 0.015 m, and 0.02 m.
3.1.1. Error measurement. In this study, the accuracy of the numerical model was evaluated using statistical indicators. The Normalized Root Mean Square Error (NRMSE) and the Mean Absolute Error (MAE) given in the following Eqs (6)-(7) were used to evaluate the models: NRMSE ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi X n i¼1 ðX i;exp À X i;sim Þ 2 n v u u t ðX exp;max À X exp;min Þ ð6Þ where n: number of data; X: hydraulic factor; exp: observed data; sim: simulated data. The lower the NRMSE and MAE values are, the better the model performance is.

Pressure head at the bottom of boxed tubes.
By using a set of piezometer tubes connected with eight gauge points along the central line of each tube, potential pressure heads were collected. Accounting for run 1, good matching was observed near the exit of the culvert, while overestimation between them was indicated near the inlet in Fig 4. Conversely, the much better agreement between numerical and experimental results was verified, as seen in Fig 5, when Q = 16.14 L/s. The finest mesh of 0.01 m gave the best result at the entrance of the intake tube, while the largest misplacement at the exit of the culvert was seen at this cell size. The coarsest one provided the largest difference between simulated and observed results ( Fig  5) and the largest values of NRMSE and MAE in both box tubes (Table 3). Moreover, the most suitable solution was provided using a resolution mesh of 0.015 m (Fig 5), and this cell size also gave the smallest values of NRMSE and MAE, indicating the most suitable grid size (Table 3). This is because the 0.015 m grid size performed better at the culvert bottom than the others; hence, it simulated more crucial pressure at the basement of the culvert. Therefore, this cell size was taken to simulate other hydraulic features and to characterize sediment scour mechanisms at the downstream channel.

PLOS ONE
Numerical study of sediment scour at meander flume outlet of boxed culvert diversion work 3.1.3. Water depth and depth-averaged velocity at the entrance of the intake. Both the measured data and numerical results showed that the free surface flow went through the twobox culvert in run 1 (Fig 6). Moreover, water depth and depth-averaged velocity at the entrance of intake measured in section 1 (Fig 2) showed the close agreement between observed and numerical solutions in all test cases, (Fig 7). In terms of the smallest input discharge, the relative error of water depth and velocity were 5.6% and 1.8%, respectively. For the largest one, 2.96% was estimated for the former and 2.38% for the later.
Consequently, these results indicate that Flow 3D is an efficient tool for the prediction of hydraulic characteristics.

Numerical simulation sediment scour in flume outlet
In this section, two modules, namely the turbulence and sediment scours of the Flow 3D package, were activated to predict scour geometry at the flume outlet. Essential computational time (t s ) was considered as maintaining steady flow, therefore obtaining equilibrium scour depth in all case studies. In subsection 3.2.1, the influence of diversion types on the relation between time t and value of d s /d s max was determined to find critical time t s , where d s is net packed sediment changed in time, and d s max is the maximum value of d s . This procedure was also described on [8] and [31] when investigating the numerical value of d s max . Then, scour behaviors after the rigid apron and at the bend segments were numerically studied.

Computational time.
Scour time is a necessary time to obtain equilibrium depth. After conducting some preliminary tests, it was observed that maximum scour depth became consistent in time when the run time reached 1000 s. Thus, we took the maximum passing discharge as a case study for this issue.
In Flow 3D, maximum scour depth (d s ) was collected every 100 seconds. The d s /d s max ratio was used to normalize scour depth, where d s max denotes the maximum scour depth. The relation between this ratio and the time corresponding to four types of wing walls was plotted in Fig 8. This figure indicated that, in the early stage, the sedimentation rate yielded by model A was greater than the others. Because, at 100 s, its d s /d s max was 75%, while values generated by models B, C, and D were 70%, 63%, and 57%, respectively. Four types of headwalls show that  all of them obtained approximately 95% of the deepest scour at 800 s (Fig 8). Consequently, 1000 s was the time required to maintain constant dimensions of scouring geometry in the packed bed of the exit channel. Therefore, in the first scouring simulation, model A was considered as the base case to verify the influence of the roughness/d 50 ratio (c s ) on scour hole geometry. The results are presented in subsection 3.2.2. engineering [19,24,32,33]. Soil information such as mean grain size (d 50 ) and roughness are important parameters affecting the dimensions of scour hole geometry. To discover the influence of these parameters on the erosion and deposition behavior in packed downstream domains, the maximum input flow, and model A of the wing wall were selected. The morphological change of downstream flow can be divided into two zones: the area after the straight section, and the bend of the meander flume. The scour behavior of each zone can be described as follows:

The influence of the roughness/d 50 ratio on the deformation of flume outlet. The necessity of reduction of scours after rigid aprons has been always a crucial issue in
• In Zone 1, the magnitude of velocity at the outlet of the culvert is much higher than that in Zone 2. The degradation of the packed sediment fume is mostly caused by the high densimetric Froude number.
• Zone 2 includes the inner and outer banks of the meander flume. Coriolis force is the main reason for deposition at the convex bank and erosion at the concave one.
Accounting for various roughness/d 50 ratios (c s ), four values of this ratio including 1.0, 1.5, 2.0, and 2.5, were used to verify the deformation of the downstream flume, in which the suggested value of c s by Flow 3D was 2.5. The development of Zone 1 caused by c s = 2.5 was much larger than that yielded by a smaller value of c s (Figs 9 and 10). Additionally, the degradation of two sides of the flume at Zone 2 was also significantly changed, although its shapes were quite similar in all cases. A strong erosion hole and a sand mound appeared on concave and convex banks together, respectively. The rate of d s max /d 50 changes in time at Zone 2 was more consistent than that at Zone 1. The greater the value of the studied ratio was, the more deeply the local scour hole at the concave bank was obtained (Figs 9 and 10).
In comparison with the empirical results taken from [5,8], and [34], the numerical solution of d s max /d 50 corresponding to four cases of c s and two values of input discharge are shown in Table 4.
• (Emami & Schleiss, 2010) [34]: where:   Numerical solution of maximum scour depth is highly dependent on the value of c s (see Table 4). Moreover, the dimension of scour hole developed strongly when the value of c s increased (see also Fig 9). Therefore, calibrating the value of c s to get the reasonable numerical solution of sediment scour must be done before simulating sediment transport.
To calibrate this parameter, three empirical equations of d s max /d 50 functioned of Densimetric Froude number (F 0 ) were selected. Table 3 indicates that the numerical results of equilibrium scour depth taken by three values of c s including 1.0, 1.5, and 2.0 are underestimated in all empirical results. While c s = 2.5 gave the suitable value with analytical solutions. Therefore, this parameter was taken to numerically assess the effect of wing wall types and its dimension on morphological change in the exit meander channel.

Effect of headwall types on erosion and deposition on channel outlet.
In the case study, four models of diversion walls at the intake outlet were investigated for their influence on morphological changes in the two zones. The scour behavior of each zone can be described as follows: • The scour hole was mostly developed due to high-speed velocity from the culvert at Zone 1.
Numerical results showed that the slope of upstream of the scour hole was much greater than that at downstream (Fig 11). In the case of no diversion wall (model B), the geometry of the erosion hole was symmetrical, while models A, C, and D generated an asymmetrical scour hole. A higher headwall prevented erosion on the concave bank, and sedimentation developed quickly on the convex bank. With the same main grain size of the rough bed, d 50 = 3.6 mm, the maximum scour depth was produced by model C, and its location occurred at the tip of the lower wing wall (Fig 12). The largest degraded area was also seen in model C and the smallest one was in model D (Fig 15). Particularly, Fig 12 delineated the deformation of sections 4 and 5. The long wing wall constructed on the concave bank of models A and C prevented this side from eroding. However, the deepest point occurred at Zone 1 was yielded by model C when the ratio d s max /d 50 obtained a maximum value of 33.9. Zone 1 had the shallowed hole created by model B when no walls were installed (Fig 12). • Moreover, the quantity of d s max /d 50 caused by all models A, B, C, and D at Zone 2 was quite similar, i.e., 29.7, 29.2, 31.7, and 28.6, respectively (Fig 13). Cross-section 5 had the same degraded form in all models. The lower erosion hole occurred on the concave bank, and the smaller one appears on the other side. The highest sand dune was observed in model A, followed by model C, while the height of the deposition mound created by models B and D were similar (Figs 12 and 13). However, Figs 13 and 14 also illustrate that at both centers of sections 4 and 5, the long headwall (case C) created the deepest erosion in both zones.
In conclusion, accounting for the sand removal volume, two short headwalls (case D) gave the smallest area affected by deposition and sedimentation. While two long ones (case C) generated the largest area and the deepest scour depth.

Effect of wing wall types on hydraulic characteristics in the meander exit channel
To explain the impact of diversion walls on erosion and deposition mechanism at two zones of the exit channel, the depth-averaged velocity distributions at cross-sections 4 and 5 were

PLOS ONE
investigated (see Figs 14 and 15). Model B (no headwalls) showed the maximum velocity appearing at the center of cross-section 4, which indicated why this case created the maximum scour depth at the center of this section. The material was removed from beneath the headwall on the right side of model C. Vector velocity fields also indicated vortex flow existing on the flow drained out of the culvert. On the other hand, all vector direction patterns in section 5 were from the concave side to the convex side, which means that fluid and soil particles may transport from the left bank to the right bank. Moreover, the streamline of depth-averaged velocity displayed in the 3D view also indicated that vortex flow appeared near the tip of the wing wall on the convex bank in all cases (Fig 15).
At the bent segment, the bed topography was strongly three-dimensional deformed (see Fig  9) [35]. The pools always formed near the outer banks, whereas the deposition fronts always formed near the inner bank. Two driving forces accounting for the gravity force and the Coriolis force are reasons of whirlpool occurred at the meander segment. Figs 16 and 17 presented the distribution of three vortex components in sections 4 and 5 in case of no headwalls (case B). At cross-section 4, all vortex components were distributed symmetrically because it was not influenced by the curved flow (Fig 16). However, in section 5, in the y-direction, clockwise and counterclockwise flows incorporated with positive and negative values of vorticity led to a whirling motion of the fluid. The higher the magnitude of the vorticity value was, the stronger whirlpool occurred (Fig 17). Hence, the soil element was detached and moved as bed-load or suspended sediment along the main flow direction and from the concave bank to the convex bank.

Conclusion
The 3D CFD model has been used widely in many complicated hydrodynamic problems. Two modules, namely turbulence and sediment scour, were selected to simulate several hydraulic characteristics of flow over a culvert coupling with a packed bed of exit flume. The numerical results of water depth, velocity and pressure head, and flow regime were compared with observed data, showing good agreement between the three operating conditions. Mesh sensitivity was analyzed with three mesh resolution sizes and indicated that a grid size of 0.015 m was adequate in simulating several sophisticated hydraulic features.
The study showed that the critical time of 1000 s was necessary to maintain a steady state for all simulation cases accounting for sediment scour phenomena. With four models of diversion walls at the outlet of the intake as well as various roughness/d 50 ratios, the influence of these scenarios on sedimentation and deposition after the straight and meander segments were discovered. The conclusions were as follows: • Two distinct zones were formed in the flume outlet. The deeper zone was formed due to high-speed velocity after the rigid apron. The erosion on the concave bank and deposition on the convex bank occurred at the meander segment due to the effects of the meandering.
• The bed deformation was highly sensitive with roughness/d 50 ratio (c s ). Particularly, the higher this ratio was, the deeper and larger the volume of bed deformation was. In comparison with empirical solutions, the numerical result of maximum scour depth yielded by value c s of 2.5 was the most reasonable.
• The influence of the four models of the wing wall on morphological change at the flume outlet was investigated. Two short headwalls caused the smallest erosion depth in Zone 1, and two long walls generated the deepest erosion hole in this area. A long wingwall on the concave bank of models A and C created the largest deformation area in this bank. Four models also had an approximately equal height of sand mounds at the convex bank. In terms of the smallest volume of bed change, two short headwalls were the best solution.

Author Contributions
Conceptualization